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Abstract 

We present a model problem for benchmarking codes that investigate 
magma migration in the Earth's interior. This system retains the es- 
sential features of more sophisticated models, yet has the advantage of 
possessing solitary wave solutions. The existence of such exact solutions 
to the nonlinear problem make it an excellent benchmark problem for 
combinations of solver algorithms. In this work, we explore a novel al- 
gorithm for computing high quality approximations of the solitary waves 
and use them to benchmark a semi-Lagrangian Crank-Nicholson scheme 
for a finite element discretization of the time dependent problem. 

1 Introduction 

Benchmark problems are of great utility for verifying and comparing nu- 
merical algorithms, and exact solutions play an important role in con- 
structing such benchmarks. Unfortunately, exact solutions may be dif- 
ficult, or impossible, to construct for nonlinear problems. In this work, 
we formulate a benchmark problem for the simplest non-linear model of 
magma migration, and explore algorithms for constructing the exact so- 
lutions and simulating the system. 

On the viscous time scale, many processes in the solid Earth, including 
mantle convection, magma migration, and crustal deformation, occur at 
such low Reynolds numbers that inert ial terms can be neglected. These 
quasi-static systems generically take the form 

dtip + V • (ipv) = sources and sinks of ip (la) 
V • cr(v; ip) = body forces on the medium (lb) 
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where ip is some material parameter, such as temperature, fluid fraction 
or chemical concentration, which influence constitutive relations such as 



permeability or rheology of the medium. Though (la) is written to suggest 
a hyperbolic nature, it may of course be parabolic, with diffusive terms 
appearing as sinks. All time dependence comes from the first equation, 



(la), which affects the velocity field through the elliptic problem (lb), 
which then advects the scalar fields. This coupling introduces its own 
set of scientific computing challenges; there is a simultaneous need for a 
robust, fast, elliptic solver and an efficient time stepping algorithm. 

1.1 A Reduced Model for Magma Migration 

An important example of a coupled hyperbolic-elliptic systems in Earth 
science is the PDE's governing the flow of a low viscosity fluid in a vis- 
cously deformable porous matrix which has been used to model the flow 
of partially molten rock (magma) in the Earth's interior. Beginning with 
the primitive equations first formulated in [15], one can, after many sim- 
plifications, arrive at the system: 

</H = ct> m V, (2a) 
[0 m - V (0 n V-)] V = -V • n e d . (2b) 

In the above equations: (i) <j> is porosity, or volume fraction of melt; (ii) V 
is the "compaction pressure", measuring viscous deformation associated 
with a compaction of the rock; (Hi) is the the unit vector in coordi- 
nate d associated with the direction of gravity. The above, dimensionless, 
Eqs. (J2| have been scaled as follows. Porosity is scaled to a reference 
value 0o ~ 0.1 — 1% and distance to the compaction length, an intrinsic 
length scale of the primitive system [l5] . The compaction length depends 
on porosity, but we scale to £o, the compaction length at the reference 
porosity. All computations here are performed with respect to the dimen- 
sionless equations, however, we will often refer to lengths in terms of the 
number of compaction lengths, multiples of So. 

Eq. |2]), as written, is more amenable to adding additional physics 
such as solid advection or melting. However, it can also be formulated as 
a single equation, together with a far field boundary condition as: 

<h + d Xd (0 n ) - V • U n V (0" m 0t)] = 0, lim 0(x, t) = 1. (3) 

|x|— too 

Derivations from primitive equations for conservation of mass, momen- 
tum, and energy are given in [7 JT8 23] . 

|3} i! 



The important feature of |3l is that it possesses solitary wave solu- 
tions, localized states which propagate in space at a fixed speed without 
changing shape. Such solutions are ideal for benchmarking as one need 
only check the distortion of the transported waveform. Indeed, in a frame 
moving with the solitary wave, the solution will appear constant. Numer- 
ical studies of the equation and its solitary waves have been performed in 
one 7,18 , two [l9], and three [33] dimensions. Recent work 20-22] has 
shown the well-posedness of ^ and the stability of its solitary waves in 
dimension d — 1. 
We demonstrate: 
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1. A novel algorithm for computing the solitary waves, which solve a 
time independent equation, based on the Cardinal Whittaker sine 
function which provide better than polynomial accuracy, 

2. A semi-Lagrangian Crank-Nicolson algorithm for a finite element 
discretization of |2|, as one algorithm for the benchmarking. 

An outline of this paper is as follows. In section[2]we review some prop- 
erties of the solitary waves. We then discuss the sine collocation method 
in section |3J"j and how it can be implemented for these equations. Then, 
in section [4] we demonstrate the the algorithm with convergence results. 
Section |5] explains the time stepping algorithm and its performance. We 
offer some remarks and comments in section [6] 



2 Solitary Wave Solutions 

The solitary wave solutions of ^ are exponentially decaying, radially 
symmetric humps in excess of <j) = 1, traveling in the Xd direction at 
a fixed speed. They are akin to the soliton solutions of the Zakharov- 
Kuznetsov equation, [34], found in plasma physics. 
Making the ansatz 



0(x,t) 0, 

the solitary waves solve the third order equation 



J2 x i + ( Xd -cty\ = </> c (r) 



(4) 



-ccj)' c + + 



+ 



1 — m 

c(d - 1) xn ( l,^.,,^, 
r 



1 — m 

-C(j)' c + Oc)' + C (0c (log C )") 

+ c(d-i)# (log </> c y 



m / 1 



m — 1 



(5) 



with the boundary conditions 



4>'M = 

lim <j> c {r) = 1 



(6a) 
(6b) 



Derivatives are taken with respect to r. 

The existence of solitary wave solutions was proven by a phase plane 
argument in dimension one. Further smoothness properties were stated 
in [21] as part of the stability analysis . The formal existence of solitary 
wave solutions in higher dimensions is an open problem, though it is 
expected. 

In one dimension, |5]) can always be integrated up twice, reducing the 
problem to quadrature. In the case case n = 3 and m = 0, after one 
integration we have 



l) + d>l-l + abWc = 



(7) 



3 



Two more integrations give the implicit expression 



r z = A + 



-2 - <j> c + 



^A^l-^AT- 



: log 



where A = (c — l)/2 is the amplitude. 



(8) 



2.1 Reformulation by Even Extension 

For higher dimensions, such integrations of |5} are not possible and we 
resort to numerical approximation. First, ^ is rewritten to make it more 
amenable to computation. Integrating from r to oo yields the integro- 
differential equation 



< 



1) + - 1) + 



1 — m 



b n M- m )") 



l) + (^-l) + c(C(log^)") 



to# 1 



TO = 1 



(9) 



The sinc-collocation method we wish to apply does not apply directly to 
problems posed on the semi-axis, (0, oo). The problem must be altered to 
live on all of (— oo, oo). Let 



4>c(x) = <t>c(\x\) 



(10) 



be the even extension of (f) c to the whole real line. Dropping the ~'s, (f) c 
solves 



i + 



±n / / 1 — m\//\ 



1 — m 



= { 



-c( ( /> c -l) + ^-l + c(^(log(/) c ) // ) 



+ c(d- 1 



J-oo 



(log( 



m / 1 



m — 1 



(ii) 



We discretize and solve (11) using the sinc-collocation method, described 
in Section l3J"1 



3 Collocation & Continuation 

Equation is a nonlinear integro- differential equation posed on an un- 
bounded domain. We employ a method, sinc-collocation, that respects 
these features. The sine spectral method is thoroughly formulated and 
explained in 12,27-29 and briefly in Appendix |A| In the sine discretiza- 
tion, the problem remains posed on R and the boundary conditions, that 
the solution vanish at =boo, are naturally incorporated. This method has 
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been applied to a variety of differential equations; see the above refer- 
ences. It has been used to compute solitary waves in the early work [l3j 
and more recently in [l4]. It was also used to study a dependent problem, 
KdV, in g]. 

Collocation insists the equation be satisfied at the nodes of the mesh. 
This is in contrast to a Galerkin formulation, which would have us dis- 
cretely orthogonalize the residual against some family of functions. Col- 
location can be interpreted as discretely solving the classical form of the 
equation, while Galerkin discretely solves the weak form. Moshen & El- 
Gamel found collocation to be superior for certain problems, [16] . 

This discretization will lead to a nonlinear system of algebraic equa- 
tions, requiring a good initial guess for the solver. Our strategy for con- 
structing such a solution is to take the d = 1 solution and then perform 
numerical continuation in dimension, up to the desired value. Construct- 
ing the d — 1 solution is also done by continuation, using an asymptotic 
approximation in the small amplitude, c ~ n, state; continuation is per- 
formed in c to its desired value. 



3.1 Sine Discretization 

Given a function u : R — »• R, u is approximated using a superposition of 
shifted and scaled sine functions: 

N N 

Cm,n(u,K)(x) = ^fcsinc ^ k ^j - ^ u k S(k,h)(x), (12) 

fc=-M k=-M 

where x k = kh for k = — M, . . . , iV are the nodes and h > 0. There are 
three parameters in this discretization, h, M, and N, determining the 
number and spacing of the lattice points. This is common to numerical 
methods posed on unbounded domains, 

A useful and important feature of this spectral method is that the 
S(k, h) functions act like discrete delta functions, 

C m ,n(u, h)(xk) = u k . (13) 

For sufficiently smooth functions, the convergence of this approximation 
is rapid both in practice and theoretically. See Theorem |A.2| in Appendix 
[A] for a statement on optimal convergence. 

Since the solution is even, we may take N — M; we write 

C M (u, h)(x) = C m ,m( u , h)(x). (14) 

We further reduce the number of free parameters down to just M, by 
slaving h to M as 

The motivation for this choice is discussed in Appendix [A] It is closely 
connected to the theory of the sine method, and the asymptotic decay 
properties of both <j) c — 1 and its Fourier transform. 
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Letting u c = <p c — 1, the solitary wave with the asymptotic state sub- 
tracted off, we formulate (11) as a nonlinear collocation problem at the 
nodes {xk}' 



+ 



-cC M (u c , h)(xk) + (Cm(u c , h)(xk) + l) n - 1 

1 — m 

c(d 



c d 2 

+ (Cm(u c , h)(x k ) + l) n -r^ (Cm(u c , /i)Ofc)) 1-m / 1 

1 — m ax z m y£ i 



r \c M (u c ,h)(x) + l) n x^ i~ [{CM{u c ,h){x)f- m ]\dx. 



-cC M (u c , h)(x k ) + (C M (u c , h)(x k ) + l) n 
d 2 

+c(C M (u c ,h)(x k ) + l) n ^ (l°g^M(w c ,/i)(xfe)) m = 1 
+c(d- 1) f°° (C M (u c ,h)(x) + l) n x^ ji A [logC M K,/i)(x)]| dx. 

" x k ^ ^ 

(16) 

From here on we suppress the subscripts c in u c and cj> c . 

Let u be the column vector associated with the sine discretization of 
u, at the collocation points {x^}, 

C M (u,h)(x k ) ^ u= (u-m u-m+1 ... u M ) T ■ (17) 

Associated with u is <fi — u + 1- 

We now define a series of matrices that operate on u. The derivatives 
of sine approximated functions are given by: 



Explicitly, 



Df k = ^SU,h)(x)\ x=Xk (18) 



w = < 


fo 














k—j 






-IT 2 
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-2(-l) fc - 









The integration matrix is 



Jk 2 7T 

where Si is the sine-integral function 



D (" 1 ) = J + ^Wj-fc)) (20) 



Si( s ) s / EEWdt. (2i) 



Jo 



Being singular, the ^ ^ operator must be treated carefully. For smooth 
even functions (V(0) = 0), it is well defined. Taking limits of the sine ap- 
proximation of an even function (u-k = Uk), 

1 7T 2 

lim -d x uoS(0,h)(x) = -wr^uo, (22) 
x— >o x on* 

lim - (d x u k S(k, h)(x) + d x u- k S(-k, h)(x)) = (23) 
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The matrix approximating is denned as: 



D 



(i) 



2(-l) fc 



I. 3h 



r j = k = 0. 



(24) 



With these matrices, the discretization of (11) is equivalent to the 
nonlinear algebraic system 



F(u) 



-cu + < 



1 + 



1 



771 



-</> n D (2) (</> 1_m 



+ 



c(d- 



^ ( - 1) (0)^ (1) 5 (1) (</> 1 - 

1 — m 

-cu + </> n - 1 + c</> n L> (2) log 

+ c(d - l)73 ( " 1) </> n J D (1) £> (1) log <f> 



m / 1 



(25) 



where 1 is a vector of size 2M + 1 with l's in all entries. Nonlinear terms 
should be interpreted as component- wise operations on the vectors. 



3.2 Initial Guesses and Numerical Continuation 



To solve (25), one needs a good initial guess of u. For dimension one, an 
excellent guess is available. Integrating |5]) reduces the equation to first 
order, which can be solved by quadrature and root finding. Sometimes it 
is even possible to obtain implicit solution, as in (JsJ . This is not possible 
for d > 1, nor is it always desirable to work out the quadrature formulas. 
Thus, for a given c and d, we proceed in two steps: 

• For d = 1, perform numerical continuation in c, from a value of 
c rsj n, up to the desired value. 

• For d > 1, apply the continuation in c to construct the solution in 
d = 1, then perform numerical condition in up to the desired 
dimension. 



3.2.1 Continuation in c 

In [32], the authors observed that in the limit of small amplitude distur- 
bances of the reference state, (3) was, to leading order, governed by the 
Korteweg - de Vries equation. Generalizing this observation in [21], let 



7= t/1 

2 (xi - ct) = 1 + 



n 
c ' 
2 



-U(j( Xl -Ct)). 



Then U solves the equation 

-U+^U 2 + diU = 0( 1 2 ), 
and small amplitude solitons, where <7< 1, are approximately 



(26a) 
(26b) 

(27) 

(28) 
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Given the desired value of c, we partition (n, c] into P points 
n < a < c 2 < . . . < cp = c 

we iteratively solve 

-cu + ct> n - 1 + j^^D^ ( < /> 1 " m - 1 



G(u;c) = 



-cu + 0"-l + I ^^ 2 » log</> 



m / 1 
m = 1 



using ir-" as the initial guess for 



G(u 



O'+i). 



;c i+ i) = 



(29) 



(30) 



(31) 



The u^ guess is given by (28). We have successfully solved with P 



O (10^), though this could likely be refined with more sophisticated con- 
tinuation algorithms. 

3.2.2 Continuation in d 

For d > 1 we solve by numerical continuation in dimension by making d 
a parameter: 



H(u;« 



-cu + < 

+ 

-cu + 



-1 + 



1 — m 
'-D { - i ^ B D (1 ^ (1, (0 1 



c ( d ~ 1) nC- 1 )^ n(!) n(!) -™ 
1 — m 



m / 1 



m = 1 



-l + c</> n L> (2) log</> 
+ c(d - l)D ( - 1) n D (1) D (1) log</> 

Given the dimension d for which we desire a solution, we partition [1, 



(32) 



into 



(33) 



1 = do < di < d2 < . . . dp — d. 
Then, assuming we have solved 

H(u w) ;d,-) = 0, 

u^ becomes the initial guess for <i/+i. P, the partition size of [1, d], need 
not be that large. P — O(lOcf) appears sufficient. As with continuation 
in c, more sophisticated continuation algorithms might improve this. 

3.3 System Size Reductions 

The even symmetry can be exploited to reduce the size of the algebraic 
system. Since u-k = Uk, we need only track Uk, k = 0, 1, . . . M. The sym- 



metry is imposed on ( 25 ) by the following manipulations on a discretized 
operator, A. Since only the last M + 1 rows are required, we only retain 
M + 1, . . . 2M + 1. Next, we add or subtract the columns 
. M onto the columns j = 2M + 1, . . . M + 2. For even/odd 
symmetry preserving operations, and we add. For even/odd 

symmetry reversing operations, ^ and J*^ , we subtract. This reduces 
the system to M + 1 points. 



for z 

-^4-2 -7 1 i — 
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c = 6, n = 3, m = 0, M = 40 +1 Pts. 
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Figure 1: Examples of computed solitary waves with different values of c, n, m, 
and M. 

4 Example Solitary Wave Computations 

We implemented our algorithm using NumPy/SciPy. The codes for com- 
puting the sinc-collocation matrices were motivated by the Matlab codes 
discussed in [31] . 

4.1 Solitary Wave Forms 

As a first example of our results, we compute a collection of solitary 
waves for different parameter values and dimensions. These wave forms 
are pictured in Figure [l] The amplitude of the wave tends to increase 
with dimension. We have not observed a choice of (c, n, m) for which this 
does not happen. This observation was previously noted in [33] for the 
n — 3, m — case. 

We can further consider this relationship between dimension and am- 
plitude by plotting "dispersion relations" between the parameter c, and 
the amplitude of the associated solitary wave. See Figure [2] The n — A. 
m — 1 case has a much more nonlinear dispersion relation than n — 3, 
m — case. 

In all cases with m — 1, we did not show any d — 3 solitary waves. 
While performing numerical continuation in d for these cases, our solver 
failed to converge at an intermediary value of d between 2 and 3. 
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Figure 2: Dispersion relations between c and solitary wave amplitude. 



4.2 Convergence 

Not only does the sinc-collocation approach work, but it also converges 
quite rapidly to the desired solution. We now present some benchmarks 
on the convergence of the amplitudes of the solitary waves in different 
dimensions, for different choices of (c, n, m). We focus on comparing the 
amplitudes (0 C (O)). Other points are more difficult to compare as the 
grids change with M. The results are given in Tables ^ [2] [3] 



5 Time Dependent Simulations 

Spectrally accurate solitary wave profiles are extremely useful as initial 
conditions for benchmarking and exploring numerical solutions of the full 
space time PDE's. For example, Figure [3] shows a solution for an off- 
center collision of two 2-D solitary waves with speeds c = 5 and c = 7 
initialized with two sinc-collocation solutions. This calculation is done in 
a moving frame translating at the mean velocity of the two waves which 
allows long runs in limited numerical domains. More specifically we solve 
a variation of the coupled hyberbolic-elliptic problem, |2]) for porosity (j) 
and compaction pressure V [10] : 

^ = <TV (34a) 
[-V • fV + <T\ V = -V • fe d (34b) 
where is the unit vector in the direction of gravity and 
D<p d<p „ A 

is the material derivative in a frame moving at speed v (which we assume 
is constant for these problems, but can vary in space and time for more 
general magma dynamics problems). 



Given any numerical method for solving (34), the sine solitary wave 
solutions provide a straightforward benchmark problem: use a high res- 
olution (M = 150) solitary wave as an initial condition and solve. A 
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Table 1: Convergence of the solitary wave amplitude for c = 4 for the 
m = problem. 



M d = 1 



20 
40 
60 
80 
100 
200 
400 
800 



d = 2 



d = 3 



1.500 21353765 
1.50000080060 



1.70417661440 
1.70608902282 



L96849289246 
1.97466312561 



1.50000000 989 1.70617 046612 1.9748 6670209 

1.500000000 23 1.70617 690685 1.97488 105982 

1.5000000000 1 1.706177 67812 1.974882 64950 

1.50000000000 1.706177828 34 1.974882937 68 

1.50000000000 1.70617782848 1.97488293789 

1.50000000000 1.70617782848 1.97488293789 



Table 2: Convergence of the solitary wave amplitude for c = 5 for the n 
m = 1 problem. 



M 


d= 1 


d = 2 


20 


14.3312283238 


22.5954364016 


40 


14.2972695906 


22.6643001828 


60 


14.2972368619 


22.6667057152 


80 


14.2972367260 


22.6668188493 


100 


14.2972367248 


22.666827544 


200 


14.2972367248 


22.6668286095 


400 


14.2972367248 


22.6668286096 



Table 3: Convergence of the solitary wave amplitude for c — 6 for the 



m = .5 problem. 



M d= 1 



20 
40 
60 
80 
100 
200 
400 



d- 



d- 



1.479 45862654 
1.47938232695 



L67961 148369 
1.68059282799 



L94941026371 
1.95217168937 



1.47938214 568 1.68062 352158 1.9522 3978609 
1.47938214411 1.68062557559 1.95224385885 



1.47938214408 
1.47938214408 



1.680625 79033 
1.68062582653 



1.952244 25280 
1.95224431473 



1.47938214408 1.68062582655 1.95224431476 
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perfect scheme would have the wave propagate at speed c with no change 
in shape; anything else is numerical error. 



5.1 Numerical Methods 



System ( 34 ) has previously been solved by finite difference and finite vol- 
ume methods with explicit time stepping and operator splitting |33| e.g.]. 
Here we describe and benchmark more recent implicit finite element codes 
with semi-Lagrangian/Crank-Nicolson time stepping for the advection 
terms. Specifically, we solve the non-linear variational problem 



F(u) = / [f n Vv • (Vp - e d ) + vf m p] dV+ f f n e d 
Jn Jan 

+ J n q(f-^f m p-9(x))dV = 



dS 



(35) 



for consistent porosity, /, and pressure, p, at time t + At. Here, u = (p, /) 
is a solution for pressure and porosity in a mixed finite element space 
V with test functions v = (v,q). For the problems shown here we use 
second order elements on triangular (2D) and tetrahedral meshes (3D) 
(i.e. V = [P2 x P2]). The semi-Lagrangian source function 

g (x*) = /(x* , t) + ^ /(x* , i)">(x* , t) , (36) 

depends on the porosity and pressure at the previous time step evaluated 
at the takeoff point of the characteristics that intersect the quadrature 
points at time t + At 25,26, e.g.]. For constant background advection 
x* = x - vAt. 



Equation (35) is non-linear with F(u) being the residual for any func- 
tion u G V. We solve F(u) = 0, using pre-conditioned Newton-Krylov 
methods implemented in hybrid FEniCS (http:/ /www. fenics.org) and[PETSc 
[3|{5] codes. FEniCS is a suite of advanced, open- source software libraries 
and applications that allows for high-level description of weak forms us- 
ing a "unified form language" (ufl) that can be translated into efficient, 
compilable C++ code using their form compiler FFC. In addition ufl has 
the capability of describing and calculating the weak form of the exact 
Jacobian (J(u) = SF/5u) by automatic functional differentiation. Given 
weak forms for both the residual and Jacobians, FEniCS provides rou- 
tines for assembly into discrete vectors and matrices, which are solved us- 
ing PETSc's non- linear equation solvers (SNES). These codes are highly 
flexible and can be used to rapidly compose a wide range of PDE based 
models and adjust solver strategies at run time. We have used these codes 
to explore a variety of physics-based block preconditioning strategies for 
iterative solutions of the magma equations. Appendix [B] provides de- 
tails of the numerical method used here. Benchmark codes are available 
through the Computational Infrastructure for Geodynamics ( CIG ) . 

5.2 Benchmark problems and results 

Table [4] gives parameters for four 2-D problems and one 3-D problem. 
For each problem we set the background advection velocity to v = — c 
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(a) 




(b) 

Figure 3: Figure showing the off-center collision of two solitary waves with wave- 
speeds c = 7 and c = 5 and material exponents n = 3, m = 0. (a) Porosity. 
(b) Compaction Pressure. Model domain is 64 x 128 compaction lengths with 
128 x 256 degrees of freedom (node spacing h = 0.5(5). Courant number is 1. 
This model is calculated in a frame moving at the mean speed of the two waves 
(V = 6) and shows the typical non-linear phase shift interaction with some 
radiation loss on collision. 
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Table 4: Parameters for solitary wave benchmarks 



n 


m c 


Amplitude 


ft (compaction lengths) 


2D 


runs: 






3 


5 


2.33407 


64 x 64 


3 


10 


5.18711 


64 x 64 


2 


1 2.5 


2.44620 


64 x 64 


2 


1 4 


11.03790 


64 x 64 


3D 


runs: 






3 


5 


2.72588 


32 x 32 x 32 



so that the wave appears stationary in the moving frame. Each model 
run is calculated on a square or cubic domain large enough so that the 
tails of the solitary wave are (f) c — 1 < 10 -7 (see Table [Z]) , with boundary 
conditions (f,p) = (1,0) on the top edge and "free-flux" (VP • n = 0) 
on the other three sides. For each wave, we consider a series of spatial 
and temporal discretizations by varying the inter-node spacing h = .25, 
0.5 and l.OJo where So is the compaction length in the constant porosity 
background. Time steps are chosen such that each wave moves a fixed 
multiple of a compaction length in a time-step i.e. cAt/5o = 0.125, 0.25, 
0.5, 1.0 and 2.0. The Courant number is cAt/h. 

Semi-Lagrangian schemes are characteristic based and do not have a 
CFL stability condition. This allows for large time steps. Nevertheless, as 
the results here show, accuracy is degraded at large time-steps. For this 
problem, the most efficient and accurate runs occur at CFL = 1. 

For these problem, two natural measures of error are distortion of 
the wave shape and disturbance of the phase speed. Given a computed 
solution for porosity /(x, £), we identify both types of errors by first min- 
imizing the functional 

£(d) = f (/(x, t) - <£ c (x + <5, t)f dx (37) 

where <p c is a high quality approximation of the solitary wave at time t. 
Though this can be done by direct minimization of the functional, it is 
more accurate to solve the nonlinear system 

<f(x)-</> c (x + «5),V</> c (x + «5))=0 (38) 

Numerical experimentation shows that this optimization/root finding 
problem reduces to a single parameter system for the vertical component 
of the displacement 8 due to the symmetries of the equation. 8 measures 
the phase shift, from which we can find the error in the speed parameter, 

8 « (c - c)t 

where c is the speed of the numerically perturbed solitary wave. Once we 
have found 8, we can directly compute the L 2 error in the approximations. 

Figure [i] shows the relative shape error \/£/||0c||2 and relative velocity 
error \c/c — 1| as a function of time for a n = 3, m = 0, c = 10 wave 
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Figure 4: Relative errors in (a) Shape and (b) Velocity for a 2-D, n = 3, 
m = 0, c = 10 wave as a function of time, grid spacing h/So and time-step 
At. Grid spacing is the distance between nodes for quadratic elements and is 
relative to the compaction length in the constant porosity background. cAt/So 
gives the number of reference compaction lengths traveled in a time step. The 
Courant number is cAt/h. For all runs, there is a rapid adjustment from the 
initial condition and then a steady evolution with little to no change in shape 
or velocity. Note, that the error depends almost entirely on the time step 
not courant number, i.e. runs with the same time step and different courant 
numbers have very similar errors. 



c=10, d=2, model=n3m0 
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Figure 5: Convergence behavior as a function of time step at t ~ 0.45 for the 
same runs shown in Fig. [4j For any given courant number, the error shows 
second order convergence and runs with integer courant number show nearly 
identical errors for the same time-step At. Runs with fractional courant num- 
ber show similar convergence but larger shape errors (and in this case smaller 
velocity errors) at the same time- step. For these runs highest accuracy per 
computational work occurs for courant number cfl = 1 runs. 
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Figure 6: Convergence behavior for all runs as a function of time- step. All 
problems show second-order convergence with At and integer courant number. 
For any given model, larger amplitude waves have larger errors and n = 2, m = 1 
waves tend to have larger errors than n = 3, m = 0, likely due to the steeper 
porosity gradients in the m = 1 waves. 3-D runs show similar errors to 2-D 
runs. Relative error for all waves is < 10 -3 for grid-spacing h = 0.25, cfl=l. 
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and a range of grid spacing and time steps. Similar to previous solitary 
wave studies [24], the initial conditions rapidly adjust to an approximate 
solitary wave in the discrete function space and then propagates with 
constant form and phase speed. The numerical methods shown here are 
second order as expected with errors that depend primarily on the time 
step 0(At 2 ) (Figures [HJj6]). Rough rules of thumb for accuracy in these 
second order codes is that grid spacing should be h < 0.25^o with time 
step given by a courant number of 1. 



6 Discussion 

One of the challenges of studying Eq. (|3| is that it is fully nonlinear. This 
carries over to the solitary wave equation |5|. As previously discussed, for 
d > 1, this cannot be reduced to a first integral. For other multidimen- 
sional equations, such as the Nonlinear Scrhodinger (NLS) equation and 
the Zakharov-Kuznetsov (ZK) equation, one can apply spectral renormal- 
ization/Petviashvilli's method [l) |ll[[T7] . However, these equations are 
semilinear, i.e., they are already of the form 

-V 2 u + \u- f(u) = 0. 

Thus, this popular approach is not immediately applicable to |5|. 

Another approach to solving for solitary wave profiles would be to use a 
shooting algorithm, though it is highly unstable, [6] |33| . Other approaches 
are comparable to our sine collocation approach: use finite differences or 
finite elements to construct a nonlinear algebraic system. The advantage 
of using sine is that it naturally incorporates the boundary conditions 
at infinity, whereas one would need to introduce an artificial boundary 
condition for these problems. 

It was necessary to extend <j) c from (0, oo) to be a function along the 
whole real axis to use sine collocation. We accomplished this by an even 
extension; however, there are other possibilities. The primary way of ap- 
plying sine methods to domains other than R is to employ a function, ip(y), 
that maps R into the domain of interest. For the semi- axis, suggested 
mappings are ip(y) — log(y) and ip(y) — log(sinh(y)), [l2j. The main dis- 
advantage to mapping, as opposed to our even extension, is that we must 
modify our scheme to accommodate the boundary condition <j)' c \ r =o — 0. 

In general, the challenge of computing a solitary wave solution in 
higher dimensions is not particular to (|3j. There is similar difficulty with 
ZK and NLS and this method is likely to prove very effective in these 
problems. 

Beyond methods for accurate computation of no n- linear waves, these 
results provide the first critical tests for any code on magma dynamics. It 



should be stressed that Eqs. (34) represent the most simplified version of 
magma dynamics that only include the contributions of non-linear perme- 
ability and rheology to porosity evolution. More general systems are re- 
quired to investigate the role of thermodynamics, chemistry, and the inter- 
action with the large scale mantle flow. However, at their core, all of these 
problems need to reproduce the non-linear waves and the benchmarks pre- 
sented here are a necessary and reasonably straightforward exercise in code 
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verification. Moreover, these larger systems, still have the general quasi- 
static structure of the basic magma equations and require accurate and ef- 
ficient multi-physics solvers for coupled hyperbolic/parabolic/elliptic sys- 
tems. The physics based, block-preconditioners demonstrated here for 
magma may be a useful approach for more general problems. 
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A Sine Approximation 

Here we briefly review sine and its properties. The texts 12,28 and the 
articles 8,27,29 provide an excellent overview. As noted, sine colloca- 
tion and Galerkin schemes have been used to solve a variety of partial 
differential equations. 

A.l Overview 

Recall the definition of sine, 



A sine approximation is only appropriate for functions satisfying cer- 
tain criteria. To define such functions we first define a strip about the 
real axis in the complex plane as 




(39) 



and for any k £ Z, h > 0, let 




(40) 



D v = {z € C | \Sz\ < v) . 



(41) 



Defining the function space: 
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Definition A.l B P (D U ) is the set of analytic functions on D v satisfying: 



\\f(t + i')\\LH-v,v) =0(\t\ a ), as t^ ±00, with a e [0,1), (42a) 
lim \\f(. + [y)\\ LP + lim \\f(.-i y )\\ LP <oo. (42b) 

y^-v y^-v 

we have the following 

Theorem A. 2 (Theorem 2.16 of [l2]) Assume f e B P {D U ), p=l,2, 
and f satisfies the decay estimate 

\f(x)\ < Cexp(-a|x|). (43) 

If h is selected such that 

h = ^pKvJ(aM) < min {ttz/, tt/^2} , (44) 

then 

- d:C M (f, h)\\ L oo < CM«" +1)/2 exp ((-V^M)) . 

This theorem justifies the sine method, and guarantees rapid convergence 
for appropriate functions. Checking that a function satisfies all the hy- 
potheses is non-trivial, and in practice it is omitted. However, it is es- 
sential to have a proper value of h to ensure convergence of algorithm; 
Theorem IA.2I states the bound on h is related to estimates of both the 
decay rate and the domain of analyticity of the function. 



A. 2 Decay Rate 

The decay rate is often easy to estimate. For the solitary waves solutions 
of |5|, the asymptotically linear equation is 



d-1 



d- 1 lf 



-c(j) c + ncj) c + c ^0' c " + 
Integrating once, 

-^(</> c -l) + ^' + ^V c = 0. (45) 
c r 

The solutions in dimensions n — 1, 2, 3 which decay as r — >• oo are 

' e -VWcr n = 1, 

1 oc { k (^/l — n/cr) n = 2, 

l e -^l-n/cr 3 
r 

ko is a Bessel function, and it decays exponentially. Continuing to assume 
that (45) governs the large r behavior, the general decay relationship is 



_ d-1 

2 (r) — 1| oc r 2 e 



(46) 



The decay rate for Theorem A. 2 is thus a — y/l — n/c. 
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A. 3 Analyticity 

The second parameter, zy, the distance into the complex plane which the 
function can be analytically continued off the real axis, is not as readily 
observed. Others have found v — tt /2 sufficient. We use v — tt/2 and can 
numerically confirm a posteriori that this is a reasonable value. 

Checking this condition is equivalent to identifying the decay rate in 
Fourier space; v satisfies 

\u(k)\ < Ke~ uM (47) 
To see this, write u{x) in terms of its Fourier Transform, 

/oo 
e txk u(k)dk. 
-oo 

Evaluating u at z = x + iy, 

/oo 
e xk e~ yk u{k)dk. 
-oo 

u(z) is defined, provided 

/oo 
e- yk \u(k)\dk < oo 
-oo 

If u is bounded by e _l/ ' fc ', then this integral is guaranteed to be finite for 

— v < y < v 

Thus, we can compute the Fourier transform and assess its decay rate to 
identify v. This was used in [30] to study the analyticity of the solutions 
to several partial differential equations. 

Using our computed solitary waves, we approximate their Fourier 
transforms as follows. First, we form the even extension by reflection. 
Then we delete the node at X-m\ the solution now sits on 2M grid points. 
This extended solution is treated as periodic on [x_m+i, xm] and its trans- 
form is computed. Figure [7] shows the resulting transforms. In these cases, 
we have resolved them to machine precision. More importantly, the solu- 
tions decay more rapidly than than e _7r//2 ' fe '. This justifies using v — tt/2. 
This value, together with a = y/l — n/c, tells us that the largest value of 
h for which we can expect convergence is 



which is ( 15 ) 



B Notes on Numerical Methods for so- 
lution of space-time PDE's 

The algorithm demonstrated here uses second-order mixed Lagrange finite 
elements (P2-P2) in space for porosity and pressure and a second-order, 
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Figure 7: Fourier transforms of the solitary waves computed by the sine 
rithm. 
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semi-Lagrangian Crank-Nicolson scheme in time. The latter approach 
discretizes the generic hyperbolic advection-reaction PDE 



£1 

Dt 



s(x, t) 



as a two level scheme via the trapezoidal rule 

/(x, t + At) = /(x* , t) + ^ [a(x, t + At) + s(x* , t)] 



(48) 



(49) 



where x* is the take-off point of the characteristic that intersects point x 
at time t + At 25,26]. We use this scheme to approximate the strong form 
of Eq. ( 34 ) and then multiply by test functions and integrate to produce 

|35i). 



the non-linear weak form of the residual at time t + At(Eq 

Semi-Lagrangian methods in finite elements can be considered a distor- 
tion and reprojection problem. For the pure advection problem,!) f / Dt = 
0, the weak- form becomes 



/ vfdx 
Jn 



I 



vf(x*)dx. 



(50) 



This is just the projection of the advected continuous function /(x*,t) 
back onto the function space. To evaluate the RHS of Eq. (50) by quadra- 
ture, /(x*) must be evaluated at the quadrature points. Therefore x* 
should be the coordinates of the take-off points of the characteristics that 
intersect the quadrature point. This is different from finite difference 
problems where the characteristics intersect the grid points. 

Given the weak form of the residual, it is straightforward to calculate 
the weak form of the Jacobian by functional differentiation. This exact 
Jacobian is assembled by FEniCS into a 2 x 2 non-symmetric block linear 
problem 

A(f) B(f,p) 
AtC(f) D(f,p) 

where Su = [Sp, 5f] T is the correction to the solution at each Newton step. 
Note if At = 0, the problem becomes linear. If we begin with an initial 
guess / = fo,p = 0, then the problem reduces to solving 



(52) 





5p 












. F f . 



(51) 



A(fo) B(f ,p) ' 




Sp 




' F p (f ) - 


M 












(34 



and 
for 



where M is the porosity mass matrix. Thus for At = 0, Sf 
5p = —A(fo)~ 1 F p (fo) which is just the discrete solution to Eq. 
the pressure given porosity fo. 

For a non-zero time-step, we solve the linear problem, Eq. (51) using 
a block-preconditioned Newton-Krylov scheme implemented in PETSc, 
using their FIELDSPLIT block preconditioned] Our preconditioner uses 



A 
AtC D 



(53) 



1 PETSc gives considerable flexibility for experimenting with a wide range of solvers and 
pre-conditioners 
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as a lower triangular block pre-conditioner with one V-cycle of algebraic 
multi-grid on the A block and 2 sweeps of SOR on the D block, and then 
GMRES on the entire Jacobian. Other choices of preconditioners and 
solvers can be passed by command line arguments. However, we find this 
particular recipe to be robust and efficient in both 2 and 3-D, converging 
quadratically in two Newton steps with total residuals ||.F(u)||2 < 10 -14 
independent of dimension, grid-size, time-step and choice of model (n, m) 
or initial condition. Using the pre-conditioner as a solver, turns the al- 
gorithm into a Picard scheme with linear convergence (and a residual 
reduction of about an order of magnitude per non-linear step). 

Given a high quality initial condition for porosity generated by the 
sinc-collocation scheme, we first project that solution onto our initial 
porosity as /o, the general time stepping algorithm is: 



Algorithm 1 SLCN-Newton-Krylov algorithm for Magma 
Require: at t = 0, step k = 0: Set initial condition fa = fo, Pk = 0, 

set At = 0, Solve for p k (k = 0) by Newton 

set At = dt 

for k = 1, 2, ... do {loop until t > tm ax } 
Set initial Guess: 
set p k =p k -i 

solve Mf k = Jq (#(x*) + At/2fj?_ lPk ) dx 
while ||jF(u)|| 2 > tol do 

Iterate preconditioned Newton-Krylov method for fk,Pk 
end while 
t^t + At 
end for 



As with all no n- linear solvers, having a good initial guess is critical to 
robustness. In the above algorithm we make a prediction for the porosity 
at the future time by solving Eq. (34i) as a projection problem with a 
lagged pressure field. M is the symmetric porosity mass matrix which can 
be solved with a few iterations of ICC preconditioned CG. 

Though the above algorithm uses a constant time step, it is straight- 
forward to implement adaptive time stepping in these two-level schemes 
for variable At. 
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